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Genetic oscillators, such as circadian clocks, are constantly perturbed by molecular noise arising 
from the small number of molecules involved in gene regulation. One of the strongest sources of 
stochasticity is the binary noise that arises from the binding of a regulatory protein to a promoter in 
the chromosomal DNA. In this study, we focus on two minimal oscillators based on activator titration 
and repressor titration to understand the key parameters that are important for oscillations and for 
overcoming binary noise. We show that the rate of unbinding from the DNA, despite traditionally 
being considered a fast parameter, needs to be slow to broaden the space of oscillatory solutions. 
The addition of multiple, independent DNA binding sites further expands the oscillatory parameter 
space for the repressor-titration oscillator and lengthens the period of both oscillators. This effect is 
a combination of increased effective delay of the unbinding kinetics due to multiple binding sites and 
increased promoter ultrasensitivity that is specific for repression. We then use stochastic simulation 
to show that multiple binding sites increase the coherence of oscillations by mitigating the binary 
noise. Slow values of DNA unbinding rate are also effective in alleviating molecular noise due to 
the increased distance from the bifurcation point. Our work demonstrates how the number of DNA 
binding sites and slow unbinding kinetics, which are often omitted in biophysical models of gene 
circuits, can have a significant impact on the temporal and stochastic dynamics of genetic oscillators. 


I. INTRODUCTION 

Genetic oscillatory networks are ubiquitous in nature 
and perform important functions. For example, the 
cell cycle oscillator regulates cell growth and division, 
whereas the circadian clock regulates the behavior of or¬ 
ganisms with respect to daily changes in light. These 
genetic oscillators are used by living systems to reliably 
coordinate various periodic internal processes with each 
other as well as with their rhythmic environment. How¬ 
ever, this presents a challenge at the cellular level because 
oscillators have to maintain proper timing (temporal co¬ 
herence of oscillation) in the presence of stochastic noise 
that arises from the small number of regulatory molecules 
in cells [1]. 

A simple mechanism to mitigate the effect of molecu- 
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FIG. 1. (a) In the Activator-Titration Circuit (ATC), the ac¬ 
tivator is constitutively produced at a constant rate and acti¬ 
vates the expression of the inhibitor, which, in turn, titrates 
the activator into inactive complex, (b) In the Repressor- 
Titration Circuit (RTC) the constitutively expressed inhibitor 
titrates the self-repressing repressor. 


lar noise would be to increase the number of molecules of 
each species [2-4]. While the number of RNAs and pro¬ 
teins made per gene can be large, most cells are funda¬ 
mentally constrained to 1-2 gene copies and are subject to 
binary noise in the first step of gene regulation (i.e., tran¬ 
scription factor binding to DNA) [5, 6]. This binary gene 
regulation noise manifests itself as a stochastic temporal 
pattern of all-or-none gene activity depending on whether 
the promoter is bound by the regulatory protein or not. 
Recent work shows that slow DNA binding/unbinding 
kinetics (also called the non-adiabatic limit) can exacer¬ 
bate the binary noise and have profound consequences 
on gene expression [7], epigenetic switching [8], and os¬ 
cillation [3, 4, 9, 10]. Faster kinetic rates and complex 
gene promoter architectures have been proposed as a way 
to suppress the effect of this binary noise. For example, 
increasing the DNA binding/unbinding rate can increase 
temporal coherence of oscillations via more precise sam¬ 
pling of the concentration of transcription factors [3, 9- 
11] or by increasing the distance from a bifurcation point 
[4, 12]. However, transcription factors often have slow 
DNA unbinding rates [13-17], which suggests that these 
mechanisms are not generally applicable. The coopera¬ 
tive binding of a transcription factor to multiple binding 
sites has also been shown to increase temporal coherence 
of oscillations [18]. However, multiple binding sites do 
not always lead to cooperativity and transcription fac¬ 
tor binding to a single DNA site may often be enough to 
effectively activate or repress transcription. 

To better understand the potential mechanisms that 
suppress the binary gene regulation noise, in particular 
the influence of slow DNA unbinding rates and multi¬ 
ple binding sites, we study an Activator-Titration Cir- 
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cuit (ATC) that has been theoretically shown to oscil¬ 
late [19]. The ATC consists of a constitutively-expressed 
activator that promotes the expression of the inhibitor, 
which then titrates the activator into an inactive het¬ 
erodimer complex (Fig. 1). Studying the ATC has two 
advantages. First, it lies at the core of animal circa¬ 
dian clocks [20] and oscillatory NF-acB signaling [21, 22] 
and has served as a model of natural genetic oscillators 
[10, 19, 23-26]. Second, the ATC generates the necessary 
nonlinearities through protein titration [27] and does not 
require cooperative binding of activator to the inhibitor 
promoter. Thus, by studying a titration-based oscilla¬ 
tor, we can better explore the kinetic effects of multiple 
binding sites on coherence independently of the effects 
that might arise from cooperativity. To obtain general 
insights that are not specific to activation, we also study 
a Repressor-Titration Circuit (RTC), which consists of 
a self-repressor and a constitutively-expressed inhibitor 
(Fig. 1). This novel titration-based oscillator is analo¬ 
gous to the ATC but uses repression instead of activation 
for the transcriptional regulation. 

We first characterize these oscillators and how they 
depend on several key parameters in Section II. We de¬ 
liberately constrain ourselves to physiological parame¬ 
ters found in a simple eukaryote S. cerevisiae, commonly 
known as budding yeast. We show that, in addition 
to slow mRNA degradation, slow DNA unbinding rates 
of transcription factors are important for providing the 
necessary delay in the negative feedback loop for oscilla¬ 
tory solutions. Thus, both the DNA unbinding rate and 
mRNA degradation rate can set the period of oscillation. 
We then demonstrate that the addition of multiple, inde¬ 
pendent binding sites has nontrivial effects on the ATC 
and the RTC. While multiple binding sites lengthen the 
period of both oscillators due to an effective increase in 
the delay of negative feedback, they dramatically increase 
the oscillatory solution space of the RTC only. This is 
because multiple, independent binding sites generate ul¬ 
trasensitivity (i.e. strong nonlinear response to changes 
in regulatory protein concentration) in repression-based 
promoters only, and thus only RTC can benefit from this 
effect. In section III, we use stochastic Gillespie simula¬ 
tions to understand the extent to which DNA unbinding 
rates and numbers of binding sites suppress the molec¬ 
ular noise in ATC and RTC oscillators. We show that 
multiple binding sites increase the temporal coherence of 
oscillations by alleviating the binary noise resulting from 
discrete gene states. We also show that slower values of 
DNA unbinding rates are best for coherent oscillations in 
simple titration-based oscillations. Last, we compare and 
contrast our results on temporal coherence with those of 
previous models of genetic oscillators in Section IV. 


on gene expression can occur transcriptionally via repres¬ 
sors [29-31] or post-transcriptionally via phosphorylation 
[32-34], degradation [32, 34, 35], or titration of activa¬ 
tors into inactive complexes by inhibitors [10, 19, 23-26]. 
The ATC is a minimal two-gene circuit that can oscillate 
through the use of protein titration both as a source of 
nonlinearity and indirect negative feedback. In the first 
phase of oscillation, high levels of free activator homo- 
dimerize, bind the promoter, and overproduce inhibitor 
until all free activator has been titrated into inactive het¬ 
erodimer. In the second phase of oscillation, the remain¬ 
ing activator will unbind from the inhibitor promoter and 
be sequestered by inhibitor, thus causing the promoter to 
return to low levels of expression of inhibitor. The levels 
of inhibitor will eventually decline to a point where free 
activator can re-accumulate and restart the cycle. 

In the RTC, protein titration is used exclusively as 
a source of nonlinearity and the negative feedback is 
directly achieved through auto-repression. In the first 
phase of oscillation, high levels of free repressor will 
homo-dimerize, bind directly to its own promoter, and re- 
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II. BIOPHYSICAL MODEL OF ATC AND RTC 

Oscillators require negative feedback with nonlinearity 
and time delays [28]. Mechanistically, negative feedback 


FIG. 2. A biophysical model for ATC (a) and RTC (b) 
with explicit transcription, translation, protein-protein and 
protein-DNA interactions. Each arrow corresponds to a reac¬ 
tion rate in Eqs. 1-3. Neither of these titration-based oscilla¬ 
tors have been built or studied by synthetic biologists. 
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press its transcription. The free repressor will be titrated 
away by the constitutively expressed inhibitor. In the 
second phase of oscillation when free repressor levels are 
low, the remaining repressor will unbind from the pro¬ 
moter, returning to high levels of transcription and the 
rapid over-production of free repressor. As we will show 
below, the indirect versus direct nature of negative feed¬ 
back in ATC and RTC is responsible for many of the 
differences between these two titration-based oscillators. 


A. ATC and RTC oscillators with a single DNA 
binding site 

Even simple genetic circuits such as the ATC and RTC 
include many reactions and parameters (Fig. 2). An ex¬ 
haustive search over all the parameter space was not fea¬ 
sible, and we decided to constrain our parameter space 
by studying synthetic gene circuits that could be built 
in budding yeast. Synthetic genetic oscillators have been 
useful tools to understand the properties of natural oscil¬ 
lators. For example, a synthetic oscillator built in bacte¬ 
ria [36] was useful in understanding entrainment capabili¬ 


ties of genetic oscillators, as well as elucidating sources of 
stochasticity that affected entrainment. Surprisingly, all 
synthetic genetic oscillators built to date have neglected 
protein titration, a common mechanism in natural os¬ 
cillators. To this end, we built a mathematical model 
of ATC and RTC oscillators using a basic leucine zip¬ 
per (bZIP) transcription factor that dimerizes and binds 
DNA, and a rationally-designed inhibitor that binds free 
bZIP into an inactive heterodimer. These synthetic com¬ 
ponents have been successfully used in yeast [37] and, 
importantly, many of the protein-protein and protein- 
DNA binding affinities of this bZIP and inhibitor pair 
are known [38, 39]; see Table I. We fixed these param¬ 
eters and scanned through a range of other biophysical 
parameters to understand which ones affect oscillation. 
Our results should help guide future experimental im¬ 
plementation of synthetic ATC and RTC oscillators in 
yeast. 

The biophysical model of our ATC and RTC circuits 
is based on chemical mass-action kinetics where the dy¬ 
namic variables are the mean concentrations of all molec¬ 
ular species. The ODEs that correspond to the reactions 
in Fig. 2 are the following: 


^ = -a[Go][X 2 ] + 9[G,] (la) 

&=a[Go][X,]-9[G,] (lb) 

^ = /3[r,] - 7 [A][/] + e 2 [XI] - S,[I] (Ic) 

® = I3[rx] - l[X][I] + e^iXI] - 27 [X ]2 + 2ei[A2] - <5p[A] (Id) 

&=^[X][I]-e,[XI]-S,[XI] (le) 

& = 7[X]2 - ei[A2] - 6p[X2] - a[Go][X2] + 9[G,] (If) 


With [rx] and [r/] described by: 

^=Po[Gt]-< 5™M (2a) 

® = pf[G,] + p,{[Gt] - [Go]) - 5^[ri] (2b) 

for the ATC, where X=A (activator), and 

^ = p/[Go] + p,{[Gt] - [Go]) - ^™[rx] (3a) 

^=Pq[Gt]-5M (3b) 

for the RTC, where X=R (repressor). The first two 
equations represent the dynamics of promoter DNA 


where [Go] and [Gi] are the mean concentrations of free 
and bound promoter, respectively. The molar concentra¬ 
tion of total DNA [Gt] = [Go]-b[Gi] = 1/{Na-V) = 1/24 
nM where Na is the Avogadro constant and V is the 
yeast cell volume; see Table 1. Here, we consider only 
a single DNA binding site, but we will later expand our 
analysis to include multiple binding sites. At any in¬ 
stant, the promoter is either free or bound. The proba¬ 
bility of free or bound promoters is equal to the ratio of 
concentrations [Go]/[Gt] or [Gi]/[Gt], respectively. The 
other equations describe the mean concentration dynam¬ 
ics of the respective molecular species such as mRNA 
(fi, rx), monomeric protein (I or X), and dimeric pro¬ 
teins {X 2 , XI), where X stands for the activator A or 
repressor R, respectively. The regulatory homodimer A 2 
associates with Go at a rate a to form Gi, which disso- 
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FIG. 3. (Color online) Parameter space of oscillatory solutions on a logarithmic scale for RTC (top) and ATC (bottom) with 
increasing DNA binding sites. The colormap shows the number of p/ values that exhibited oscillations for each combination 
of /, Sm, and 8. (a) and (b), single DNA binding site for RTC and ATC, (c) and (d) three independent binding sites for RTC 
and ATC, (e) and (f), three synergistic binding sites for RTC and ATC, where activation/repression strength is when more 
than one activator/repressor is bound. 


dates at the rate 0. The r/ and rx are the inhibitor and 
activator/repressor mRNAs. For the ATC, the activa¬ 
tor niRNA [rx] is transcribed const it utively at the rate 
Po, where as the inhibitor mRNA is transcribed at rates 
Pf and pb from free and bound DNA, respectively (Eqs. 
2a,b). In contrast, for the RTC, the repressor mRNA is 
transcribed at rates pf and pb (Eqs. 3a,b) while inhibitor 
mRNA r/ is constitutively transcribed at the rate po- We 
assume that all mRNA species are degraded at the same 
rate 5m and translated into proteins with the same rate 
j3. The activator/repressor X protein dimerizes into ac¬ 
tive homodimer X 2 and forms inactive heterodimer XI 
with the inhibitor protein I at the same rate 7 . The ho¬ 
modimer and heterodimer dissociation rates are ei and 
€ 2 , respectively. We assume that all protein species are 
stable and diluted by cell growth at rate 5p. 


B. DNA unbinding kinetics influence oscillation 

Our parameters were restricted to physiological val¬ 
ues from yeast (see Table I). Most parameters were kept 
fixed, but we varied four key parameters. The first pa¬ 
rameter was the mRNA production rate (p/) of free, un¬ 
bound promoter because a desired expression level can 
easily be selected from existing promoter libraries [40]. 
Second, we varied the activation/repression strength (/), 
which is the ratio of the larger p divided by the smaller 
p. Thus, / = pb/pf for the ATC and / = Pf/pb for the 


RTC. The ratio / can be tuned by appropriate choice of 
activation or repression domains fused to our bZIP tran¬ 
scription factor [41-43]. The third parameter was the 
mRNA degradation rate (Sm), which is known to set the 
time scale of the ATC oscillator [19]. Last, we varied the 
DNA unbinding rate (6) because it is our point of focus 
and this parameter can vary between different transcrip¬ 
tion factors. The DNA dissociation constant (Kd) fixes 
the DNA binding rate a = ^; see Appendix for details. 

We divided the physiological range of each variable pa¬ 
rameter into 30 values (on a logarithmic scale) and evalu¬ 
ated the long-term dynamics of a total of (30)^ parameter 
sets per circuit. We solved the ODEs over time for each 
set oi {pf, f, 6m, 0)- A solution was classified as oscil¬ 
latory if the trough of activator or repressor homodimer 
concentration was below the Kd of DNA-binding and if 
the peak was above 2Kd; see Appendix for justification. 
We noticed that p f had the smallest effect on the number 
of oscillatory solutions and, thus, we plot the marginal 
frequency distribution of oscillatory solutions over f, 5m, 
and 0 in Fig. 3. We see that strong activators (large / for 
the ATC), stable mRNAs (small 5m), and slow DNA un¬ 
binding rates (small 0) generally favor oscillation. The 
last two parameters dictate the timescale of the delay 
in the negative feedback loop. Increased delay supports 
oscillation and, thus, the largest number of oscillatory so¬ 
lutions occur at the smallest 0 and 6m for both RTC (Fig. 
3a) and ATC (Fig. 3b). The parameter space of stable 
oscillations is larger in ATC relative to RTC for a single 
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FIG. 4. The DNA unbinding rate 9 sets the period of the 
oscillations for RTC (a) and ATC (b) at slow unbinding rates. 
The mean period of oscillatory solutions for a given 6 is shown 
(solid black line) with the shaded area representing the range 
of periods. 


binding site because of the additional step and delay in 
the negative feedback loop: negative feedback through 
the activator in the ATC is indirect (i.e. activator regu¬ 
lates the expression of inhibitor, which then inhibits its 
activity), where as the self-repressor in the RTC is direct 
(i.e. repressor regulates its own expression). 

The period of oscillation r should be set by the 
timescale of the slowest parameters in the delay. The 
negative feedback in our circuits is dominated by DNA 
unbinding rate 6 and mRNA degradation rate 5m [19]. 
This relationship can be seen in Figure 4 where the DNA 
unbinding rate sets the oscillation period at low 6. An 
increase in 9 leads to mRNA degradation rate {5m) be¬ 
coming the slower timescale at which point r becomes 
flatter and less dependent on 6. Eventually, a bifurca¬ 
tion occurs at a critical, maximum value of 0max which 
leads to loss of the stable limit cycle. A similar relation¬ 
ship exists for the mRNA degradation rate 5m] see Figure 
S2. 


Pb Pb 



Same Transcription Rate 


FIG. 5. Transitions between promoter states for multiple 
DNA binding sites (n=3). Gi denotes the set of promoter 
states with i out of total n binding sites occupied by activa¬ 
tor (X=A) or repressor (X=R) dimers. There are n — i ways of 
switching from Gi to [Gi+i] via the binding of X 2 , and there 
are i ways of switching from Gi to Gi_i by the unbinding of 
X 2 . Our model conservatively assumes that the binding of 
X 2 does not affect the binding or unbinding of the next tran¬ 
scription factor to an adjacent site (no cooperativity). We 
also assume that the transcription rate is equal to pb for Gi, 
G 2 , ... Gn promoter states. 


that if the occupancy of any binding site by a transcrip¬ 
tion factor activates or represses transcription, then the 
effective unbinding rate ( 0 „) from a state of saturated 
DNA binding to the unbound DNA state (Go, where the 
transcription rate changes) should decrease with the in¬ 
creasing number of binding sites (u). We can show that 
On = 0/Hn, where is the n-th harmonic number (see 

Supplement [44]). 

The addition of multiple DNA binding sites to our 
model will modify Eqs. (la-b. If) by increasing the num¬ 
ber of promoter states that can be bound by X 2 ] see Fig¬ 
ure 5 and Supplement [44]. For three binding sites (n=3), 
our new Eqs (la-b) are: 

^ = -3a-[Go][X2]+0[Gi] 

^ = 3a • [Go][X 2 ] - (0 + 2 a • [A 2 ])[Gi] + 29 ■ [G 2 ] 

^ = 2 a • [Gi][X 2 ] -{29 +a- [A 2 ])[G 2 ] + 30 • [G 3 ] 

^=a-[G2][X2]-30-[G3] (4) 

where the total concentration of DNA [Gt] = [Go] -b 
[Gi] -b ... -b [Gn] is fixed to 1/{Na ■ V) = 1/24 nM. For 
three binding sites (n=3), the term —a[Go][Ar 2 ] -b 0[Gi] 
in Eq. (If) is replaced with: 


C. Multiple DNA binding sites affect ATC and 
RTC oscillators differently 

This role of DNA unbinding rate in generating de¬ 
lays led us to hypothesize that multiple DNA binding 
sites should increase the parameter space of oscillations 
and lengthen the period of the oscillators. We reasoned 


3 3 

-^(3-*)a-[G.][A2]+^*-0[G,] (5) 

2—0 2—0 

Gi denotes the promoter states with i out of total n 
binding sites occupied by activator (X=A) or repressor 
(X=R) dimers. The factors in front of each term repre¬ 
sents the amount of degeneracy of each state, i.e [Gi] has 











6 


i bound sites, thus i ways of switching to There¬ 

fore, we have the term i- 6 [Gi]. At the same time, [Gi] has 
n — i vacant sites, so it has n — i ways of switching to the 
state [Gi+i]. Thus, we have the term (n — i)a ■ [Gi][A 2 ]. 

The addition of two more independent DNA binding 
sites dramatically increased the oscillatory space of the 
RTC (Fig. 3c), while slightly decreasing the oscillatory 
space of the ATC (Fig. 3d). These opposite results arise 
from a compound effect. First, two extra binding sites 
decreased the effective unbinding rate for the promoter 
to be fully vacated by half {O/H 3 « 9/2). This decrease 
in effective 9 increased the delay and resulted in some im¬ 
provement in oscillations for both RTC and ATC. This 
effect is best observed in the increased period of both 
oscillators (Fig. 4). The second, more dominant ef¬ 
fect is the fundamental difference in how the promoter 
sensitivity changes with multiple, independent binding 
sites. It is well established that nonlinear promoter re¬ 
sponses facilitate oscillation [28]. We use the logarithmic 
sensitivity (S) to quantify the nonlinearity in the pro¬ 
moter response, where S = dlogP/dlog[A 2 ] [45]. P is 
the synthesis rate of the promoter and [X 2 ] is the acti¬ 
vator/repressor homodimer concentration that regulates 
the promoter. As shown previously [45, 46], an increase 
in the number of independent repressive binding sites will 
increase the magnitude of S and create an ultrasensi¬ 
tive promoter response, (i.e. [S'] > 1, see Supplement 

[44] ). However, increasing the number of independent 
activating binding sites cannot generate an ultrasensi¬ 
tive promoter response ([S'] < 1); see Discussion in [45]. 
In fact, the logarithmic sensitivity for activation actu¬ 
ally decreased with the number of binding sites at our 
physiological concentrations (see Supplement [44]). This 
difference is the reason why the RTC and ATC oscillators 
exhibited fundamental differences to increased number of 
binding sites. Our work shows that synthetic repression- 
based oscillators are preferable designs because the RTC 
gets an effective boost in promoter ultrasensitivity simply 
by adding multiple, independent binding sites. 

We also tested whether synergistic repression or activa¬ 
tion might change our results. Synergistic activation or 
repression occurs when the states that have more than 
one binding site occupied (i.e. G 2 and G 3 ) are acti¬ 
vated/repressed /^-fold instead of /-fold because they 
can interact with RNA polymerase at several interfaces 

[45] . Although this synergy increased the activation or 
repression strength, it did not significantly change the 
oscillatory parameter space (Fig. 3e,f). 


III. STOCHASTIC SIMULATIONS 

Deterministic simulations were useful for understand¬ 
ing how DNA unbinding rate and the number of bind¬ 
ing sites affect the phase space and period of oscillation. 
However, they cannot provide insights into the loss of 
temporal coherence that arises from stochastic gene ex¬ 
pression. To this end, we used the Gillespie algorithm 



Time (min) Deiay (min) 


FIG. 6. Sample stochastic trajectories for one and three bind¬ 
ing sites for RTC (a) and ATC (c), and their autocorrelation 
functions (b,d). The variable parameter values {6, 5m, Pf, f) 
were fixed to (0.02 min ^, 0.0159 min 0.8928 min 3.63) 
for the RTC and (0.02 min“^, 0.0186 min“^, 0.1781 min“^, 
30) for the ATC. We chose parameters that produced oscilla¬ 
tion over the largest range of DNA unbinding rates. The rest 
of the parameters are given in Table 1. 


[47] to simulate stochastic chemical reaction kinetics and 
investigate how DNA binding/ unbinding dynamics and 
the addition of binding sites affect the temporal coher¬ 
ence of ATC and RTC oscillators. 

For each ATC and RTC, we quantified temporal co¬ 
herence by calculating the autocorrelation function of 
mRNA transcripts levels (repressor mRNA for the RTC 
and inhibitor mRNA for the ATC); see Fig. 6 . In the ab¬ 
sence of noise, an undamped oscillatory signal will have 
an undamped, periodic autocorrelation function. The 
presence of noise will stochastically perturb period and 
phase, such that the autocorrelation now exhibits damp¬ 
ening or loss of temporal coherence. We quantihed the 
loss of coherence by measuring the rate of exponential 
decay (e“*/'^“) of the envelope of a periodic (cos(27rt/r)) 
autocorrelation function (see Appendix for details). Sim¬ 
ilar to other studies [4, 18], our metric for temporal co¬ 
herence is the normalized autocorrelation function decay 
rate, which is the ratio of timescales Tq/t. A larger ra¬ 
tio indicates better temporal coherence. We varied the 
DNA unbinding rate (9) and number of binding sites (n) 
to understand the role of each feature in resisting molec¬ 
ular noise. 


A. DNA Unbinding Rate 

Our results show that ATC and RTC oscillators with 
smaller DNA unbinding rates exhibit better temporal co- 
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Sites Sites 

FIG. 7. The normalized autocorrelation function decay rate 
for the RTC (a) and ATC (b) for varying 6 and number of 
binding sites. All parameters, except 0, are the same as in 
Fig. 6. Boxed, outlined regions are parameters past the 
bifurcation point (0max) where deterministic oscillations are 
unsustainable and damped, yet exhibit stochastic excitable 
oscillations. 



herence (Fig. 7). Lower 9 increases the temporal co¬ 
herence of the oscillations because of the increased dis¬ 
tance of the dynamical system from the bifurcation point 
(^max); see Figure 8. Eventually there is another bifur¬ 
cation at small dmin, but these unbinding rates are un- 
physiological and do not affect our conclusions regard¬ 
ing biophysical ATC and RTC oscillators. Strikingly, 
some 9 > dmax, which do not show deterministic oscil¬ 
lation, exhibit oscillation in the presence of noise. This 
phenomenon is consistent with coherence resonance [48] 
which has been observed in other excitable, genetic cir¬ 
cuits near oscillatory bifurcation points [12, 24]. 


B. Multiple DNA Binding Sites 

Increasing the number of binding sites (n) also in¬ 
creased the temporal coherence of ATC and RTC oscil¬ 
lators over all DNA unbinding rates (Fig. 7). To bet¬ 
ter understand this result, we must consider the effect 
of stochastic binding and unbinding of regulators on the 
variance of gene expression. In the phase of changing ac¬ 
tivator/repressor concentrations, the binding sites start 
being occupied or vacated. Each additional binding site 
introduces an additional DNA binding state. Because 
we treat the expression level of all bound DNA states as 
equivalent (pt,), the spontaneous binding and unbinding 
events that occur between states that have at least one 
binding site occupied have no effect on transcription; see 
Eig. 5. These “protected” states act as a buffering mech¬ 
anism to mitigate the effects of binary noise on temporal 
coherence. 


IV. DISCUSSION 

We analyzed the properties of two titration based ge¬ 
netic oscillators, the activator-titration circuit (ATC) 
and the repressor-titration circuit (RTC). The focus of 
our study was to understand how the number of DNA 


(a) 




0 (min ’) 


FIG. 8. Bifurcation diagram of the RTC (a) and ATC (b) 
oscillators as a function of DNA unbinding rate (6). All pa¬ 
rameters, except 9, are the same as in Fig. 6. There are two 
bifurcation points (0max, 6^min) and the amplitude of mRNA 
oscillation is shown by the upper and lower branches. Phys¬ 
iological values of 6 are to the right of the dashed vertical 
line. 


binding sites and slow unbinding kinetics in promoters 
mitigate or exacerbate the binary gene regulation noise. 
First, we showed that multiple DNA binding sites and 
slow unbinding kinetics were important for providing the 
necessary delay in the negative feedback loop for oscilla¬ 
tory solutions. The role of slow DNA binding/unbinding 
in providing delay for oscillations is consistent with prior 
work on a small negative feedback oscillator [9]. Sec¬ 
ond, we used stochastic simulation to show that slower 
DNA unbinding rates exhibited better temporal coher¬ 
ence, a result which appears at odds with previous work 
on circadian clocks and NF-kB oscillators [3, 10, 12] and 
which is more in line with the results obtained for a small 
negative feedback oscillator model [9]. This previous 
work showed that slower DNA unbinding kinetics neg¬ 
atively affected the temporal coherence for two reasons. 
First, slow DNA unbinding increased the stochasticity of 
gene expression due to imprecise concentration sampling, 
which decreased the temporal coherence of oscillations 
[3, 10]. Second, slower DNA unbinding (0) pushed the 
dynamical system towards dmin bifurcation point, which 
made it less robust to noise [12]. These results are differ¬ 
ent from ours because the delays in the circadian clock 
and NF-kB models rely on slow intermediate steps (e.g. 
phosphorylation and/or nuclear transport) in the nega¬ 
tive feedback loop. Unlike our titration-based oscillators 
in Fig. 8, these models do not have ^max and still oscil¬ 
late at infinitely fast unbinding rates where the promoter 
dynamics are in steady-state. 

We observed the opposite effect for our titration-based 








oscillators because physiological 6 overlaps the ^max bi¬ 
furcation point for ATC/RTC. Thus, lowering 0 always 
increases the robustness in ATC/RTC because the dy¬ 
namical system is moving away from dmax and deeper into 
oscillatory parameter space. This phenomenon likely ex¬ 
plains the similar results presented in [9]. The influence 
of DNA unbinding rate on temporal coherence depends 
on the structure of the underlying bifurcation diagram of 
each oscillator as a function of 9. Changes in topology, 
mechanism, and parameters can change the bifurcation 
diagram and, thus, the influence of DNA unbinding rate 
on temporal coherence of oscillation may also change. 

Last, we demonstrate that multiple independent bind¬ 
ing sites consistently increased the temporal coherence 
of oscillations by alleviating the binary noise result¬ 
ing from binary gene states. Our results agree with 
previous work, which showed that multiple, coopera¬ 
tive DNA binding sites increased the coherence of cir¬ 
cadian clocks [18]. However, in contrast to our results, 
the temporal coherence of circadian clocks peaked at 3 
binding sites and then decreased with additional sites. 
The difference likely arises from our slower DNA bind¬ 
ing/unbinding rates, where the ATC and RTC oscilla¬ 
tors spend significant time in protected states that are 
buffered against molecular noise. In contrast, the circa¬ 
dian clock model spends very little time in the interme¬ 
diate protected states between fully free or fully bound 
promoters and, therefore, the increased coherence is only 
due to cooperativity [18]. The idea of buffering to re¬ 
duce noise in gene circuits has been discussed in the con¬ 
text of decoy binding sites [49]. However, this requires 
fast DNA binding/unbinding, where as buffering through 
promoter states requires slow DNA binding/unbinding. 
We note that increased temporal coherence due to pro¬ 
tected states is a stochastic effect because the addition 
of binding sites consistently increased the coherence of 
ATC oscillators, despite occasionally pushing it past the 
bifurcation point at dmax (Fig. 6b). 


Appendix A: Parameter Values 

To constrain the physiological parameters of our mod¬ 
els, we used data from large-scale studies of the yeast 
transcriptome and proteome; see Table I. These data pro¬ 
vide typical ranges and values for our parameters. First, 
we converted numbers of molecules into nanomolar (nM) 
concentrations using the cell volume V = 40 fL for hap¬ 
loid yeast. For the ATC, we assumed that the basal 
mRNA transcription would be low. Thus, p/ for the ATC 
was constrained to values from the bottom 5th percentile 
to the median of all mRNA synthesis rates [50]. Simi¬ 
larly, pf for the RTC was constrained to values from the 
median to top 95th percentile. In the case of the ATC, 
the constraint pf < po < Pb ensured that the inhibitor 
can completely titrate the constitutively-expressed acti¬ 
vator when the inhibitor is maximally produced at pb, but 
not when it is expressed at the basal rate p/. Similarly, 


for the RTC, the constraint Pb < Po < Pf ensures that 
constitutively-expressed inhibitor can completely titrate 
the repressor when the repressor is produced at the re¬ 
pressed rate p&, but not p/. We set po = y/PfPb to satisfy 
both conditions. mRNA degradation rate ranged from 
the bottom 5th percentile to top 95th percentile values 
for all genes [50]. To obtain a rough approximation for 
the translation rates, we assumed a constitutive gene ex¬ 
pression model for all genes: 

^=P-^rn[r] (Al) 

&=(3[r]-6,[P] (A2) 

At steady state, /3 = Protein concentrations 

and degradation rates were taken from [51, 52]. We cal¬ 
culated P for all genes and used the median value in our 
model. We also assumed that our activators, repressors, 
and inhibitors are not actively degraded and are diluted 
by growth. Thus, 5p = ln{2)/T, where T = 90 mins is 
the duration of the yeast cell cycle. The proteins in our 
models were based on a mammalian transcription fac¬ 
tor basic leucine zipper (bZIP) protein C/EBPa and its 
dominant-negative inhibitor (3HF) [37]. We used previ¬ 
ously measured rates for protein-protein interaction ki¬ 
netics [37]. Since we did not know the DNA unbind¬ 
ing rate for C/EBPa, we considered the range for the 
known DNA unbinding rates for other bZIP proteins 
[13-16]. The thermodynamic dissociation constant {Kj) 
of C/EBPa to its specihc DNA binding site is known 
[39]. We set the DNA association rate to a = 9/Kd. 
Einally, we varied the activation/repression strength / 
from 1 to 30, to consider both strong and weak activa¬ 
tion / repression. 

Appendix B: Methods 

We scanned the parameter space for oscillations by 
running simulations on MATLAB (Mathworks) using 
odel5s for 2000 min and recording the minima and max¬ 
ima of the activator/repressor homodimer during the last 
1000 min. We imposed the following restrictions: 1) the 
last minimum should be below Kd so that DNA-binding 
is not saturated, 2) the last maximum should be greater 
than 2 Kd so that the change in transcription is noticeably 
altered. While this restriction slightly underestimates the 
number oscillatory solutions, it ensures that a synthetic 
version of these circuits would produce detectable oscilla¬ 
tions. We verihed that our definition gave similar results 
to a less stringent criterion for oscillation. 

We used the direct Gillespie method to perform the 
stochastic simulations [53]. We ran the simulations for 
10 ® min and recorded the concentration of the regu¬ 
lated mRNA (inhibitor for the ATC and the repres¬ 
sor for the RTC). We then normalized the concentra¬ 
tion such that the average would be zero and evaluated 
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TABLE I. Parameter values 


Parameter 

Min 

Max 

Reference 

6l(min-i) 

0.0188 

34.5 

[13-16] 

ATC p/(min“^) 

0.0509 

0.1781 

[50] 

RTC p/(min“^) 

0.1781 

0.8928 

[50] 

(5m(min~^) 

0.0159 

0.1516 

[50] 

/ 

1 

30 


ATC pi,(min“^) 

/ • Pf 


RTC p(,(min“^) 

Pf/f 


a(nM“^ min“^) 

0/3.344 nM 

[39] 

po(min"^) 

VPfPb 



14.1 

[50-52] 

(5p(min“^) 

0.0077 


7 (nM min“^) 

0.6 

[37] 

ei(min“^) 

6 

[37] 

e2(min“^) 

0.024 

[37] 

U(fL) 

40 

[27] 


the autocorrelation function. We then fit the function 
C{t) = ■ cos(27rt/r) to the first 1500 min of the 

autocorrelation function to measure the decay constant 
To and period r. The ratio tq/t describes how rapidly 
the envelope of autocorrelation function decays per oscil¬ 
lation period. 
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I. PROMOTER WITH MULTIPLE BINDING SITES 


We first derive the chemical kinetic equations for 2 binding sites and then generalize to 
n binding sites [1]. Similar to previous work [2], we calculate the promoter sensitivity of 
activators and repressors in the limit of thermodynamic equilibrium. 


A. Promoter with 2 equal binding sites 

We hrst begin with an explicit derivation that includes each binding state combination: 


—LJ = —a ■ [Goo] [^ 2 ] + ^ ■ [Gio] — OL ■ [Goo][W 2 ] + Q ■ [Gqi] (1) 

—= a ■ [Goo] [^ 2 ] ~ ^ ■ [Goi] — a ■ [Goi][W 2 ] + 0 ■ [Gn] (2) 

—= a ■ [Goo] [^ 2 ] ~ ■ [Gio] — oi ■ [Gio][W 2 ] + 6 ■ [Gn] (3) 

^ ^ = Oi ■ [Gio][W 2 ] — 6 ■ [Gii] + a ■ [Goi][W 2 ] — 6 ■ [Gn] (4) 


where [G^] = [Goo] + [Goi] + [Gio] + [Gn] and [X 2 ] is the homodimer transcription factor 
concentration. We rewrite these equations in a compact form: 


— — 2q; ■ [Go][W2] + 0[Gi] (5) 

—LJ = 2q; ■ [Go][W2] — 0[Gi] — q;[Gi][X2] + 29[G2\ (6) 

^ = a.|GJ[X2]-2«|GJ (7) 

where [Go] = [Goo], [Gi] = [Goi] + [Gio], and [G2] = [Gn]. This simplihcation is valid 
when all binding sites have equivalent kinetics and the regulator X 2 binds independently to 
each site. 


i. Promoter synthesis rate at equilibrium 
At equilibrium, Eqs. (5-7) are equal to 0. Thus, 


2 



( 8 ) 

( 9 ) 

( 10 ) 


2|G„]|X2] = Kj ■ |Gi] 

{Kd + [^2])[Gi] = 2[Go] [^ 2 ] + ‘^Kd ■ [G 2 ] 

[Gi] ■ [X 2 ] = 2Kd ■ [G 2 ] 

where Kd = 0/a. Solving these equations with constraint [Gt] = [Go] + [Gi] + [G 2 ]: 


[Go] 

[Gi] 

[G2] 


[Gt] 

[Gt] 

[Gt] 


1 


(1 

+ 

X2] \2 

Kd ! 

2 

. ( 

[^ 2 ]n 


\ 

Kd > 

(i 

+ 

IK2] \2 
Kd ^ 

1 

r PGI'i 2 

1 

^ R 

'd ' 

a 

+ 

[^2] \2 
Kd ^ 


( 11 ) 

( 12 ) 

(13) 


The promoter synthesis rate P 2 where any bound state has the same gene expression [pb) 
and the unbound state has gene expression (pf) is given by: 


^ 2 (^ 2 ) 

where X 2 = 


[Gt] 


Pb [2 ■ X 2 + (X2)^] + pf 

( 1 +X 2)2 


[Gt] 


Pb + 


jPf - Pb) 

(1 + X2)2 


(14) 


B. Promoter with n equal binding sites 

By generalizing our previous derivation for an arbitrary number of binding sites, we get 
n + 1 differential equations, where n is the total number of binding sites and the index i 
spans 1 < i < n — 1. 


d[Go] 

dt 


= —na 


[Go] [^2] + 0 [Gi] 


(15) 


— (n — i + 1)q; ■ [Gj_i][X 2 ] — {i ■ 6 + {n — i)a ■ [X 2 ])[Gi] + (f + 1)6 ■ [Gj+i] (16) 
^ = a ■ [Gn-i][X 2 ] - nO ■ [Gn] (17) 
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Gi denotes the promoter states with i out of total n binding sites occupied by activator 
(X=A) or repressor (X=R) dimers. The factors in front of each term represents the amount 
of degeneracy of each state, i.e [Gj] has i bound sites, thus i ways of switching to [Gj_i]. 
Therefore, we have the term i ■ d[Gj\. At the same time, [G*] has n — i vacant sites, so it has 
n — i ways of switching to the state [Gj+i]. Thus, we have the term (n — i)a ■ [X 2 ][Gj]. Note 
that our model assumes that the binding of a transcription factor does not affect the binding 
or unbinding of the next transcription factor to an adjacent site (no cooperativity). We also 
assume that transcription is fully activated/repressed for any bound state of the promoter, 
i.e. transcription is equal to pb for Gi, G 2 , ..., G„. Although not shown, the ODE for the 
rate of change in dimer concentration Eq. (If) in the main text is also modihed to reflect 
binding and unbinding from each state G*. 


1. Promoter synthesis rate at equilibrium 

Setting Equations (15-17) equal to zero and solving with the constraint [Gt] = [Go] -l- 
[Gi] -|- ... -T [G„]: 


[Go] — [Gt] 


[G.] = [Gt] 

[Gtv] = [Gt] 


1 


(i + 7fr 


V! ('lAdv 

i\{N-i)\ \ Ki ) 

(l+l*)" 


([^\N 
^ Ka > 


(i + @ 


N 


The promoter synthesis rate is thus given by: 


(18) 


(19) 


( 20 ) 


P„(V2) 


[Gr] ■ 


Pb + 


(P/ - Pb) 

{I+X 2 Y 


( 21 ) 


where X 2 = 

^ Kd 


4 



C. Logarithmic sensitivity analysis 


We compute the logarithmic sensitivity of the promoter: S 


dlogP 

dlogX2 


X 2 dP . 
P(X2) dX2 ■ 


SiX^) 


_ n(pt - Pf)X2 _ 

(1+X2).(p/ + Pi.|(1 + A'2)”-1]) 


( 22 ) 


1. Activation 

We consider the simplihed case of ideal activation, where p/ = 0. Thus, 

S' is a monotonically decreasing function of X 2 with boundaries S'(O) = 1 and S'(cx)) = 
0. Thus, activation is never ultrasensitive (i.e. IS"! < 1 for any value of X 2 , n, or pf,). 
Additionally, IS"! is also a monotonically decreasing function of n for all X 2 > 0; see Figure 
SI. The boundaries are at |S'(0)| = (^i^x 2 'yin{i+X 2 ) ~ 0- This explains why the 

ATC oscillator phase space decreases with the addition of multiple DNA binding sites. 


2. Repression 

We consider the case of ideal repression where = 0. Thus, 


^(X2) 


-nX2 

(I + X2) 


(24) 


S' is a monotonically decreasing function of R with boundaries S'(O) = 0 and S'(oo) = —n. 
Promoter repression is ultrasensitive (i.e. [S'! > 1) for X 2 > and n > 1; see Figure SI. 

I S'I is a monotonically increasing function of n for all X 2 > 0. The boundaries are at 
|S'(0)| = 0 and |S'(oo)| = 00 . Thus, the RTC oscillator phase space increases with the 
addition of multiple DNA binding site. 


II. KINETICS OF COMPLETE UNBINDING FROM MULTIPLE SITES 

The probability of unbinding for a single binding site is given by a Poisson distribution: 
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pit) = 6e 


(25) 


-et 


where 6 is the unbinding rate. The average unbinding time for a single site is given by 
< r >= t ■pit)dt and is equal to Consider the case where all n binding sites are hlled 
and where the concentration of regulator abruptly changes to 0 (such that no rebinding 
occurs). We calculate the probability distribution of time for complete unbinding from all 
n binding sites. The probability that the slowest unbinding event is r and the other n — 1 
events are less than r is: 


'P(r|n, 6) 
'P(r|n, 9) 



When we normalize the distribution (Jg°° P(r|n, = 1): 


(26) 

(27) 


■p(r|n,9) = n [1 - ' ■ fe-"' 


(28) 


The average < r > can be calculated by r ■ P(r)(ir: 


< r >= 


9 9r, 


(29) 


where Hn = l + ^ + | + ... + ^ (Harmonic number) and 9n is the effective DNA unbinding 
rate. 
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III. SUPPORTING FIGURES 


3 . 

2 . 

1 . 

0 . 

FIG. SI. Multiple binding sites affect activation- and repression-based promoters differently. The 
absolute value of logarithmic sensitivity S is plotted as a function of number of binding sites N, 
with X 2 = 2. The points of 1 and 3 binding sites considered in this study are marked with circles. 
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FIG. S2. The mRNA degradation rate 5m sets the period of the oscillations for RTG (a) and ATG 
(b) at slow mRNA degradation rates. The mean period of oscillatory solutions for a given 5m is 
shown (solid black line) with the shaded area representing the range of periods. 
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